0.1 Temperaturas mínimas en España

Objetivo de aprendizaje

Este caso práctico muestra como leer y graficar datos espaciales en R. Para ello, vamos a trabajar con los datos de temperatura mínima registradas en España por las estaciones metereológicas de la AEMET.

Tarea 1: Abrimos RStudio

El presente análisis se va a realizar empleando RStudio, por lo que empezaremos abriendo el programa y creando un nuevo script de R en File>New File>R script.

Tarea 2: Importamos y describimos los datos objeto de estudio

El conjunto de datos proporcionado (tempmin.csv) contiene el nivel de temperatura del aire en España entre el 6 y el 10 de Enero de 20211. Estos datos han sido descargados usando la librería climaemet (Pizarro, Hernangómez, & Fernández-Avilés, 2021) y han sido posteriormente tratados para su uso en esta práctica.

El primer paso consiste en importar la base de datos de temperatura mínima. El archivo está en formato csv, por lo cual es un fichero de texto plano. Podemos usar varias funciones para realizar la importación, en este caso vamos a emplear paquetes del tidyverse para realizar todo el tratamiento de datos:


library(readr)

tmin <- read_csv("data/tempmin.csv")

Q1: ¿Qué información tiene tmin?

Table 1: Detalle del objeto tmin
fecha indicativo tmin longitud latitud
2021-01-06 4358X -4.7 -5.880556 38.95556
2021-01-06 4220X -7.0 -4.616389 39.08861
2021-01-06 6106X 4.7 -4.748333 37.02944
2021-01-06 9698U -6.8 0.865278 42.20528
2021-01-06 4410X -3.4 -6.385556 38.91583
2021-01-06 1331A 1.0 -7.031389 43.52472

Podemos observar que la tabla generada contiene 5 columnas distintas:

Tarea 3. Convertir tmin a un objeto de la clase espacial geoR

Para convertir un objeto a geodata (el formato requerido por geoR), proporcionaremos una tabla con las coordenadas y los valores a incluir en cada coordenada. En este ejemplo, vamos a emplear sólamente los datos correspondientes al 8 de enero.

library(dplyr)
library(geoR)

tmin_geoR <- tmin %>%
  filter(fecha == "2021-01-08") %>%
  # Seleccionamos las columnas de interés
  dplyr::select(longitud, latitud, tmin) %>%
  # Y creamos el objeto geodata
  as.geodata(
    coords.col = 1:2,
    data.col = 3
  )


summary(tmin_geoR)
#> Number of data points: 211 
#> 
#> Coordinates summary
#>      longitud  latitud
#> min -9.291389 35.27639
#> max  4.215556 43.78611
#> 
#> Distance summary
#>         min         max 
#>  0.01024389 13.85144264 
#> 
#> Data summary
#>        Min.     1st Qu.      Median        Mean     3rd Qu.        Max. 
#> -14.9000000  -4.6000000  -0.5000000  -0.6293839   3.5000000  13.6000000
plot(tmin_geoR)
Objetos en geoR

Figure 1: Objetos en geoR

DIEGO, y yo por qué se que España es esto 4326…

Está explicado en la sección de CRS, es el código del CRS que usan los GPS, el habitual para coordenadas longitud/latitud

Tarea 4. Convertir tmin a un objeto de la clase espacial sf

En esta tarea, convertiremos los datos de tmin en un objeto espacial sf, es decir, datos espaciales de tipo vector.

Los datos de tmin contienen coordenadas geográficas longitud/latitud, asi que como vimos en la sección Sistema de Referencia de Coordenadas (CRS) el CRS a emplear ha de ser un CRS geográfico. Usaremos el código EPSG 4326, que corresponde a coordenadas geográficas y suele ser el habitual en este tipo de situaciones.

library(sf)

tmin_sf <- st_as_sf(tmin,
  coords = c("longitud", "latitud"),
  crs = 4326
)

tmin_sf
#> Simple feature collection with 1066 features and 3 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -9.291389 ymin: 35.27639 xmax: 4.215556 ymax: 43.78611
#> Geodetic CRS:  WGS 84
#> # A tibble: 1,066 x 4
#>    fecha      indicativo  tmin             geometry
#>  * <date>     <chr>      <dbl>          <POINT [°]>
#>  1 2021-01-06 4358X       -4.7 (-5.880556 38.95556)
#>  2 2021-01-06 4220X       -7   (-4.616389 39.08861)
#>  3 2021-01-06 6106X        4.7 (-4.748333 37.02944)
#>  4 2021-01-06 9698U       -6.8  (0.865278 42.20528)
#>  5 2021-01-06 4410X       -3.4 (-6.385556 38.91583)
#>  6 2021-01-06 1331A        1   (-7.031389 43.52472)
#>  7 2021-01-06 1690A       -0.1 (-7.859722 42.32528)
#>  8 2021-01-06 8489X       -8   (-0.255833 40.43333)
#>  9 2021-01-06 8025         2    (-0.494167 38.3725)
#> 10 2021-01-06 9784P      -10       (0.224722 42.63)
#> # ... with 1,056 more rows

Tarea 5: Dibujemos las estaciones de monitoreo de la temperaria mínima en un mapa de España. Ámbito espacial.

Vamos además a incluir una capa de las comunidades autónomas de España. Para ello usaremos un paquete API que nos proporciona esta información en formato sf:

library(mapSpain)
# sf object
esp <- esp_get_ccaa() %>%
  # No vamos a usar Canarias en este análisis
  filter(ine.ccaa.name != "Canarias")


plot(esp$geometry) # Dibujamo el mapa de España menos las Islas Canarias
Mapa de España (Sin Canarias)

Figure 2: Mapa de España (Sin Canarias)

Q2: ¿Tengo el Sistema de referencia de coordenadas (CRS) de las estaciones de monitoreo en la misma proyección que el contorno de España?

Como se comentó en la sección correspondiente, cuando se emplean datos geográficos provenientes de varias fuentes, es necesario asegurarse de que ambos objetos están usando el mismo CRS. Vamos a comprobarlo:

st_crs(tmin_sf) == st_crs(esp)
#> [1] FALSE

Vemos que no lo están, por lo que vamos a proyectar las coordenadas a un CRS común. En este caso usaremos el CRS de referencia de tmin_sf:

esp2 <- st_transform(esp, st_crs(tmin_sf))

st_crs(tmin_sf) == st_crs(esp2)
#> [1] TRUE

Dibujamos las estaciones de monitoreo con el contorno de España. Vamos a usar el paquete ggplot2 como referencia, sin embargo existen varios paquetes especializados en mapas temáticos, como pueden ser tmapo mapsf.

library(ggplot2)

ggplot(esp2) +
  # Para graficar objetos sf debemos usar geom_sf()
  geom_sf() +
  geom_sf(data = tmin_sf) +
  theme_light() +
  labs(
    title = "Estaciones de monitoreo AEMET en  España",
    subtitle = "excluyendo las Islas Canarias"
  ) +
  theme(
    plot.title = element_text(
      size = 12,
      face = "bold"
    ),
    plot.subtitle = element_text(
      size = 8,
      face = "italic"
    )
  )
Estaciones de AEMET en la Península Ibérica

Figure 3: Estaciones de AEMET en la Península Ibérica

Q3. Mis datos y el contorno de España están en el mismo CRS, pero ¿es el adecuado?

DIEGO, POR QUÉ TRANFSORMAMOS, PARA PASAR A UTM E INTERPRETAR EN METROS?

Explicado más adelante, para hacer el grid en metros cuadrados. Si no, se haría en grados cuadrados (¿?) porque long/lat no es medida de superficie

Tarea 6: Representamos la variable temperatura mínima tmin para el día 8 de enero de 2021.

En la siguiente tarea, seleccionaremos los datos correspondientes al 8 de enero de 2021 y crearemos un mapa temático en el que representaremos los valores de temperatura mínima registrados en cada estación mediante un código de colores.

tmin_8enero <- tmin_sf %>%
  # seleccionamos el día y la variable
  filter(fecha == "2021-01-08")


plot(tmin_8enero["tmin"],
  main = "Temperatura mínima (8-enero-2021)",
  pch = 8
)
Mapa de puntos con temperatura mínima

Figure 4: Mapa de puntos con temperatura mínima

Podemos utilizar el ámbito espacial, el contorno de España para graficar y contar la historia de la Filomena un poco mejor.

Mapa completo con temperatura mínima

Figure 5: Mapa completo con temperatura mínima

Q4. El mapa ha quedado muy claro. Vemos como los datos nos cuentan la historia de Filomena en aquellos sitios donde se tomaron mediciones, pero ¿podríamos tener un mapa de interpolación para tener una estimación de la temperatura mínima en las partes donde la AEMET no tiene estación de monitoreo?

Tal y como se avanzó en teoría, parece lógico pensar que aquellos puntos que estén cerca tendrán valores similares así que tomemos ventaja de la dependencia espacial y utilicemos un método determinista, como la Distancia Inversa Ponderada, comúnmente conocido por su acrónimo inglés IDW (Inverse distance weighted), el cual es uno de los métodos más simples para llevar para llevar a cabo una interpolación espacial.

En este tipo de análisis, es crucial que el CRS sea el apropiado. En este caso, ya definimos el CRS como un CRS geográfico (es decir, usando coordenadas de longitud y latitud). Sin embargo, para el ejercicio de interpolación es más adecuado usar un CRS local (que provoque pocas deformaciones en la proyección de España) y en alguna unidad de distancia, como metros (ya vimos que en los CRS geográficos las unidades son grados).

Si usamos el paquete crsuggest podemos observar los CRS sugeridos:

library(crsuggest)

sugiere <- suggest_crs(tmin_8enero, units = "m", limit = 5)

# Usamos la sugerencia del paquete
crs_sugerido <- st_crs(sugiere[1, ]$crs_proj4)

esp3 <- st_transform(esp2, crs_sugerido)
tmin_8enero3 <- st_transform(tmin_8enero, crs_sugerido)

A continuación, para realizar la interpolación, necesitamos generar una malla que representará las celdas de las que queremos obtener el valor interpolado.

Dado que hemos proyectado nuestros datos a un CRS cuya unidad son los metros, podemos definir el tamaño de cada celda en metros cuadrados. En este caso vamos a usar celdas de 49 kms cuadrados (7 x 7 kms):

malla_sf <- st_make_grid(
  esp3,
  cellsize = 7000
)

Graficamos la superficie para ver exactamente qué hemos construido:


ggplot(esp3) +
  geom_sf() +
  geom_sf(
    data = malla_sf,
    size = 0.1,
    col = "red", alpha = 1,
    fill = NA
  ) +
  geom_sf(
    data = tmin_8enero3,
    aes(fill = "AEMET Stations"), size = 5, shape = 21,
    color = "blue"
  ) +
  scale_fill_manual(values = adjustcolor("blue", alpha.f = 0.2)) +
  theme_void() +
  theme(legend.position = "bottom") +
  labs(
    title = "Cuadrícula espacial para interpolar",
    fill = ""
  )
Malla de puntos para interpolación

Figure 6: Malla de puntos para interpolación

Se puede observar claramente cada una de las celdas que hemos creado. La interpolación asignará un valor a cada uno de ellas.

A continuación podemos llevar a cabo la interpolación usando el paquete gstat. Además, en lugar de celdas (polígonos) es necesario usar puntos en la interpolación. Calcularemos por tanto un punto representativo de cada celda, el centroide, que es el punto resultante de realizar la media arimética de las coordenadas de loss puntos que componen los lados de cada celda

# Calculamos centroide
malla_sf_cent <- st_centroid(malla_sf, of_largest_polygon = TRUE)

library(gstat)
tmin_idw <- idw(
  # Indicamos la variable que queremos interpolar
  tmin ~ 1,
  # Indicamos el conjunto de datos donde está la variable
  tmin_8enero3,
  # Indicamos la malla de destino, en sf
  newdata = malla_sf_cent,
  idp = 2.0 # Especifica la potencia de la IDW
)
#> [inverse distance weighted interpolation]
head(tmin_idw)
#> Simple feature collection with 6 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 145790.9 ymin: 68957.31 xmax: 180790.9 ymax: 68957.31
#> CRS:           +proj=lcc +lat_1=40 +lat_0=40 +lon_0=0 +k_0=0.9988085293 +x_0=600000 +y_0=600000 +a=6378298.3 +rf=294.73 +pm=madrid +units=m +no_defs
#>   var1.pred var1.var                  geometry
#> 1  2.513414       NA POINT (145790.9 68957.31)
#> 2  2.572912       NA POINT (152790.9 68957.31)
#> 3  2.633632       NA POINT (159790.9 68957.31)
#> 4  2.695589       NA POINT (166790.9 68957.31)
#> 5  2.758802       NA POINT (173790.9 68957.31)
#> 6  2.823285       NA POINT (180790.9 68957.31)

Representamos la interpolación con un mapa y mapa de contorno muy utilizado para representar datos espaciales. Para ello, vamos a usar el paquete raster convirtiendo nuestro objeto interpolado.

# Convertimos de sf a SpatiaPixels
# Esto funciona porque nuestros puntos sf están espaciados regularmente

tmin_pixels <- tmin_idw %>%
  as("Spatial") %>%
  as("SpatialPixels")


library(raster)
# Creamos un raster de nuestros pixels
rast_esp <- raster(tmin_pixels)

# Transferimos valores del objeto sf al raster
rast_esp2 <- rasterize(
  tmin_idw,
  rast_esp,
  field = "var1.pred",
  fun = mean
)

# Además, podemos recortar el raster a la forma de España

rast_esp_mask <- mask(rast_esp2, esp3)

plot(rast_esp_mask, col = colores)
contour(rast_esp2, add = TRUE)
Mapa raster con lineas de nivel

Figure 7: Mapa raster con lineas de nivel

Podemos realizar el mismo mapa usando ggplot2 y la función geom_contour_filled:


# Creo una tabla para geom contour
coordenadas <- st_coordinates(tmin_idw)
valor <- tmin_idw$var1.pred

idw_df <- data.frame(
  # Necesitamos redondear las coordenadas
  latitud = round(coordenadas[, 2], 6),
  longitud = round(coordenadas[, 1], 6),
  tmin = valor
)

ggplot() +
  geom_contour_filled(
    data = idw_df,
    aes(x = longitud, y = latitud, z = tmin),
    na.rm = TRUE,
    breaks = cortes
  ) +
  # Reajustamos la escala de colores
  scale_fill_manual(values = colores) +
  # CCAA
  geom_sf(data = esp3, fill = NA) +
  theme_minimal() +
  theme(axis.title = element_blank()) +
  labs(
    fill = "Temp. (º)",
    title = "Temperatura mínima interpolada",
    subtitle = "8 de Enero 2021",
    caption = "Datos: AEMET"
  )
Mapa con ggplot2 y lineas de nivel

Figure 8: Mapa con ggplot2 y lineas de nivel

Representamos en 3D el mapa anterior, muy utilizado también en datos espaciales

DIEGO, se útiliza mucho en espacial pero no se si este caso es muy ilustrativo….r

Por mi lo quitamos

persp(rast_esp_mask, theta = 2000, d = 0.1)

Pizarro, M., Hernangómez, D., & Fernández-Avilés, G. (2021). Climaemet: Climate AEMET tools. https://doi.org/10.5281/zenodo.5205573

  1. Las fechas seleccionadas coinciden con el periodo en el que la tormenta Filomena tuvo su auge en la Península Ibérica.↩︎